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ABSTRACT 

We present ^^Barolo, a new code that derives rotation curves of galaxies from emission¬ 
line observations. This software fits 3D tilted-ring models to spectroscopic data-cubes and 
can be used with a variety of observations: from Hi and molecular lines to optical/IR recom¬ 
bination lines. We describe the structure of the main algorithm and show that it performs 
much better than the standard 2D approach on velocity fields. A number of successful 
applications, from high to very low spatial resolution data are presented and discussed. 
^^Barolo can recover the true rotation curve and estimate the intrinsic velocity dispersion 
even in barely resolved galaxies (~2 resolution elements) provided that the signal to noise of 
the data is larger that 2-3. It can also be run automatically thanks to its source-detection and 
first-estimate modules, which make it suitable for the analysis of large 3D datasets. These 
features make ^^Barolo a uniquely useful tool to derive reliable kinematics for both local and 
high-redshift galaxies from a variety of different instruments including the new-generation 
IFUs, ALMA and the SKA pathfinders. 
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1 INTRODUCTION 


The study of the gas kinematics is one of the most powerful tools to 
infer the distribution of dark and luminous matter and understand 
the dynamical structure of disc galaxies. In this context, spectro¬ 
scopy plays a primary role since it can naturally trace the motion of 
matter thanks to the Doppler shift of the emission lines. In partic¬ 
ular, observations of the neutral hydrogen (Hi) emission-line (e.g. 
[Rogstad & Shostak|1971[|Allen, Goss, & van Woerden|1973] l allow 

to analyze the kinematics of galaxies and derive rotation curves out 
to very large radii. The total mass and the gravitational potential can 
be inferred by decomposing the rotation curve into the contribution 
given by the different matter components in a galaxy, i.e. the gas, 
the stars and the dark matter (DM) (e.g. |van Albada et al.|1985| l. 
Observed rotation curves broadly have two different shapes (e.g. 
[Casertano & van Gorkom|19M) : large spiral galaxies show steeply 
rising rotation curves that quickly become and stay flat in the ex¬ 
ternal regions. Dwarf galaxies have slowly rising rotation curves 
that may not reach a flat part. The former are believed to be domin¬ 
ated by baryons in the inner regions and by DM in the outer regions. 
The latter are thought to be DM-dominated at all radii and they have 
been largely used to constrain the inner slope of the density profile 
of DM halos (e.g. K uzio de Naray, McGaugh & de Blok|2008[|van| 
[Eymeren et al.|2009y 


* ^^Barolo is available for download at http: //editeodoro. github. 
io/Bbarolo 


Rotation curves are also fundamental to study the most im¬ 
portant scaling relations for spiral galaxies: the “Tully-Fisher Re- 
lation” (TFR, |Tully & Fisher||1977| ) and the “Baryonic TFR” 
( [McGaugh et al.|2000| l, correMing the rotation velocity (hence the 
dynamical mass) to the stellar luminosity and to the total baryonic 
mass, respectively. The TFR is a severe testing ground for galaxy 
formation: current observations show a slope L (or Mbar) oc V® with 
cr ~ 4 (e.g. |Verheijen]|2001[ |McGaugh||2012j > and these relations 
strongly constrain models in the standard ACDM picture (e.g. |Mo,| 
[Mao & White|1998] [Bullock et al.|200l] l. Finally, gas kinematics 
can also be used to infer the distribution of angular momentum in 
galaxies and study how it evolves with time (e.g. [Romanowsky &[ 
[Fall[2012[ |. 


All these fields of study require a precise determination of the 
rotation curves, from the innermost to the outermost parts of galaxy 
discs. The most common approach to describe the kinematics of 
spiral galaxies make use of the so-called “tilted-ring model” (e.g. 
[Rogstad, Lockhart, & Wright|1974[ |. Such a model is based on the 
assumption that the emitting material is confined to a thin disc and 
that the kinematics is dominated by the rotational motion, so that 
each ring has a constant circular velocity Vrot(^). depending only on 
the distance R from the centre. This approximation is very good for 
local disc galaxies without bars as their orbits are nearly circular, 
with a typical axis ratio bla > 0.9 (e.g.[Franx & de Zeeuw|1992[ 

|Schoenmakers, Franx, & de Zeeuw|1997| l. In a tilted-ring model, 

the disc is therefore broken down into a number of concentric rings 
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with different radii, inclinations, position angles and rotation velo¬ 
cities. 

The standard approach to link the tilted-ring model to observa¬ 
tions is through the fitting of 2D velocity fields. The observed velo¬ 
cities along the rings can be written as a finite number of harmonic 
terms (e.g. [Franx, van Gorkom, & de Zeeuwjp^^ . At the first 
order and neglecting the streaming motions, this expansion leads 
to the expression V\os{R) - Fsys + VrotW cos ^ sin/, where Vsys is 
the systemic velocity and 6 is the azimuthal angle, measured in the 
plane of the galaxy {6 -Q for major axis), related to the inclination i 
and the position angle 0, measured on the plane of the sky (see Fig. 
1^. It is then straightforward to fit such a function to a velocity field 
by using any of the numerous non-linear least-squares fitting tech¬ 
niques (e.g. [Press et al.|20()7| Chap. 15). There are several available 
algorithms that perform this fit to velocity fields. The most used 
is Rotcur ( [van Albada et al.|1985[ |Begeman|1987|), whic h can be 
found in data analysis packages like AIPS ( |Fomalont|1981| ), GIPSY 
( |van der Hulst et ar]|1992| ) and NEMO ( |Teuben||1995| ). Rotcur 
simply fits the above function by using a Levenberg-Marquardt 
solver { [Levenbergl 1 944[|Marquardt| 1 963| ). An optional term for the 
expansion/contraction motions through the disc is also available in 
Rotcur. In the last ten years, several new codes have refined this ba¬ 
sic 2D approach by experimenting different fitting techniques and 
by extending the harmonic expansion to higher-order terms to take 
into account asymmetries and non-circular motions. These 2D fit¬ 
ting algorithms include Reswr i (|Schoenmakers|199^ , Ringfit PU 
[mon_e^jT][200^, Kinemetry ([IG^noyiu_e^^]J200^ and DiskFit 
( [Spekkens & Sellwood|2007l|Sellwood & Sanchez|2010| ). 

All the above 2D algorithms are fast from a computational 
point of view and do return good kinematic models and reliable 
rotation curves when applied to high-resolution velocity fields. A 
number of drawbacks however exist. First of all, a velocity field is 
not the full data set, but it is itself derived from a data-cube, requir¬ 
ing the intermediate step of extracting a characteristic velocity from 
the line profile at each spatial pixel. Most used approaches include 
intensity-weighted mean along the velocity axis moment, Rog-| 
|stad & Shostak|197r|l, fi tting single ( |Begeman|1987[[Swaters| 999| l 
or multiple ( |Oh et al.||2008| ) Gaussian functions to the profiles, 
and/or including a skewness term (e.g. hs Gauss-Hermite polyno¬ 
mial, [van der Marel & Franx|199^ . Unfortunately velocity fields 
derived with different methods can significantly deviate from each 
other, especially for galaxies with asymmetric profiles. Moreover, 
the derivation of an unambiguous velocity is not possible when the 
line of sight with respect to the observer intersects the disc twice 
or more, like in nearly edge-on systems or in the presence of thick 
discs and outer flares. 

The most severe problem in deriving the kinematics from the 
velocity field is however the beam smearing (e.g. [Bosma l[T978l 
|Begeman|1987^ . The finite beam size of a telescope causes the line 
emission to be smeared on the adjacent regions. The effect is that 
the gradients in the velocity fields tend to become flatter while the 
line-profiles in each spatial pixel become broader. In other words, 
part of the rotation velocity is turned into line broadening that can 
be erroneously interpreted as gas velocity dispersion, producing a 
well-known degeneracy between these two quantities (see Section 
|3.3| ). The typical effect of the beam smearing is that the derived 
rotation curves will rise very slowly in the inner regions with po¬ 
tentially dramatic consequences for the interpretation (e.g. [Lein:] 
[Fraternali & Sancisi|2010| |. The effect becomes more and more pro¬ 
nounced as the spatial resolution of the observations decreases and 
the inclination angle of the galaxy increases. 

The natural solution for almost all the above-mentioned issues 


is a 3D modeling of the data-cubes. A 3D approach gives larger op¬ 
portunities for modeling galaxies and is effectively not affected by 
the beam smearing, since the instrumental effect is introduced in 
the model through a convolution step (see Section |2^. The main 
drawback is the computational slowness. Unlike the 2D tilted-ring 
model, an analytic form for the fitting function in 3D does not exist 
and the model is instead constructed by a Monte-Carlo extraction 
(see Section [TT) . Thus, the fit must be performed with algorithms 
that do not require the knowledge of any partial derivative (e.g. 
[Press et al.|2007| Chap. 10). Such techniques are known to be com¬ 
putationally expensive and they may converge to a local minimum 
of the function. In addition, the minimization is not just performed 
on a single map, but on the whole data-cube, which consists of n 
maps, where n is the number of channels. Finally, a larger number 
of parameters than in 2D is needed to describe the model. 

The visual comparison between the data-cube and an artifi¬ 
cial model cube has been used as an additional step to improve the 
results of the 2D tilted-ring model and potentially correct for the 
beam smearing, especially in the investigation of the shape of dark 
matter profiles in the inner region of galaxies (e.g. [Swaters|1999[ 
[Gentile et al.|2004| l. A pioneering attempt to directly use a 3D ap¬ 
proach was made by [Corbelli & Schneider[ \\991\ in a study on 
the outer warp of M33. A currently available algorithm that can 
directly fit a 3D tilted-ring model to data-cubes is TiRiFiC ( [J6zsa[ 
[et al.[2007] l. TiRiFiC has reached a considerable degree of develop¬ 
ment and sophistication and it has been successfully used to study 
the kinematics of nearby galaxies with peculiar features, like strong 
warps and extra-planar gas ( [Jozsa et ST][2009[ [Zschaechner et al.j 
[2012[ [Gentile et al.[[2013|. A 3D approa ch is also used in a re¬ 
cent tool (GalPak3D, [Bouche et al.|2015] l targeted to high-redshift 
galaxies. 

In this paper, we present a new software to automatically fit 
3D tilted-ring models to emission-line data-cubes. The name of the 
code is ^^Barolo, which stands for “3D-Based Analysis of Rotat¬ 
ing Objects from Line Observations”. This code works with 3D 
FITS images having two spatial dimensions and one spectral di¬ 
mension (velocity, frequency or wavelength). ^^Barolo builds a 
number of models in the form of artificial 3D observations and 
compares them with the input cube, finding the set of geometrical 
and kinematical parameters that better describes the data. Our pur¬ 
pose is to provide an easy-to-use algorithm that might be applicable 
to a wide range of emission-line observations, from radio-Hi data of 
local galaxies to sub-mm and optical/IR lines of galaxies up to high 
redshift, from high to very low spatial resolution. Unlike TiRiFiC, 
which has been mainly developed to study local galaxies with a 
detailed description of peculiarities, such as warps, spiral arms and 
lopsidedness, ^^Barolo is designed to work on low-resolution data, 
where the kinematic information is largely biased by the size of the 
beam. Therefore, we kept the disc model as simple as possible and 
we restrain the number of parameters to the bare minimum. 

This paper is organized as follows. In Section we illustrate 
^^Barolo’s main algorithm, going through the main fitting steps 
and describing their features. In Sectionwe show some applica¬ 
tions and tests with Hi data. We start with the modeling of a well- 
known local galaxy at high spatial resolution and then move to 
lower resolution data. We discuss the robustness and the limits of 
^^Barolo and we show the differences between a traditional 2D 
and our 3D approach in data largely afflicted by beam smearing. In 
SectionHwe test the accuracy of the code by using mock galaxies. 
Section [^recaps and discusses possible developments and future 
applications of ^^Barolo. 


© 2015 RAS, MNRAS 000, [I]-?? 























































^^Barolo; a new 3D algorithm to derive rotation curves of galaxies 



Figure 1. Flowchart of the ^^Barolo main algorithm. For each ring R, the 
algorithm builds a 3D model, convolves it with the observational beam/PSF 
and compares it with the data. If no convergence has been achieved, 
^^Barolo updates the parameters and starts over. When the algorithm con¬ 
verges to the minimum, it moves to the next ring. The optional normaliza¬ 
tion step takes place after the convolution step. 

2 THE MAIN ALGORITHM 

The main feature of ^^Barolo is that it simulates data-cube ob¬ 
servations, starting from the model of a rotating gaseous disc, and 
compares them with real data. The disc is made up by a number 
of concentric rings with non-zero thickness. The emission from gas 
in each ring is generated in a 6D domain (three dimensions for the 
spatial location and three for the components of the velocity) and 
these rings are then projected into the 3D data-cube (i.e. two spa¬ 
tial and one spectral dimensions). The comparison with the data is 
perfomed ring by ring. At each step, if the model is good enough 
(see below), the algorithm moves to the following ring, otherwise 
it updates the disc parameters until the best match between model 
and observations is found. 

The minimization between model and data is performed by 
using the multidimensional downhill simplex solver, also known 
as the Nelder-Mead method ( |Nelder & Mead||1965| ) for the min¬ 
imization of non-analytic functions. The user supplies a number of 
initial guesses for each model ring from which the function to be 
passed to the minimization algorithm is built. Input parameters can 
be either globally or ring-by-ring defined. If no guess is given, the 
algorithm will autonomously estimate the initial parameters for the 
fit (see Section [O] ). The function to minimize is an indicator of 
how the model and the data differ from each other and its construc¬ 
tion occurs through four main steps: 

(i) Disc model. The disc model is built by a Monte-Carlo re¬ 
production of the gas distribution both in the space and velocity 
domain. This function derives from the Galmod routine ( |Sicking| 
\\991) implemented in GIPSY. 


3 

(ii) Convolution. The model is degraded to the same spatial 
resolution of the data via the convolution with a 2D Gaussian rep¬ 
resenting the observational Point Spread Function (PSF). The spec¬ 
tral broadening is instead taken into account in the model construc¬ 
tion (see Section [TT] ) 

(iii) Normalization. The model is normalized to the 0th mo¬ 
ment map of the observations pixel-by-pixel or azimuthally. This 
step can be skipped and the density profile supplied by the user or 
possibly fitted. 

(iv) Residuals. Comparison between the model and the data 
pixel-by-pixel. The sum of the residuals is returned to the minim¬ 
izing algorithm. 

A schematic flowchart of the main fitting algorithm is shown 
in Fig.[^ In the following sections, we describe the most important 
steps and the main features of the algorithm. 

2.1 Disc model 

The artificial gaseous disc is constructed with a 3D tilted-ring 
model by using a stochastic function that randomly populates the 
space with emitting gas clouds, for which line profiles are built. 
Each ring of radius R and width W, is described by the following 
geometrical and kinematical parameters: 

- Spatial coordinates of the centre (aq, yo). 

- Systemic velocity Vsys- 

- Inclination angle i w.r.t. the observer (90° for edge-on). 

- Position angle 0 of the major axis on the receding half of the 
galaxy, taken anticlockwise from the North direction on the sky. 

- Rotational velocity Vrot- 

- Velocity dispersion cTgas. 

- Face-on gas column density E. 

- Scale-height of the gas layer zo- 

All these quantities are allowed to vary from ring to ring. The 
first six parameters are the same required by 2D fitting algorithms 
like Rotcur. The geometry of the tilted-ring model is shown in 
Fig. Each ring is filled with gas clouds whose spatial position 
is given in cylindrical coordinates by a radius Rc (with R-WI2 < Rc 
< R+W/2), an azimuthal angle 6c (0 ^ 6c < In) and a height Zc 
above the plane of the disc. Radius and azimuth are randomly and 
uniformly chosen, the height is selected as a random deviate from 
a given vertical distribution of the gas density (available functions 
are Gaussian, sech^, exponential, Lorentzian and box layer). The 
position of the clouds is then rotated and projected onto the plane 
of the sky with a given orientation with respect to the observer, 
according to the position angle and inclination at that radius. 

Once the positions of the clouds are determined, the observed 
velocities along the line of sight are calculated as a combination 
of systemic, rotational and random motions. The velocity profile at 
each location is built around the average velocity by dividing the 
clouds into a number of sub-clouds with velocities distributed as a 
Gaussian with dispersion cr^ = being cTgas the intrinsic 

gas dispersion and crinstr the instrumental broadening. These velo¬ 
cities are then discretized and the contribution of the sub-clouds 
is recorded in a model cube with the same sizes of the data-cube. 
^^Barolo uses by default crinstr = tVch/ V2 In 2, where Wch is the 
channel width of the data-cube. This is usually a good assumption 
in Hi data-cubes where Hanning smoothing has been applied ( |Ver-| 
|heijen|19W| l. Otherwise, the user can supply an own value for the 
spectral resolution. 
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2.2 Convolution 

The model has a nominal spatial resolution of a pixel and it needs 
to be smoothed to the same spatial resolution of the data. This re¬ 
quires to perform a spatial convolution by the observational PSF, or 
beam, for each spectral channel. We approximate the PSF as a two- 
dimensional Gaussian function, which is an adequate choice for 
radio, millimeter/submillimeter and optical/IR observations. The 
Gaussian function is defined by three parameters that character¬ 
ize the elliptical cross-sectional shape of the kernel: the full-width 
half-maximum (FWHM) of both the major and minor axes (a and 
b, respectively), and the position angle of the major axis {xj/), meas¬ 
ured anticlockwise from the vertical direction. The kernel of the 
two-dimensional Gaussian is defined by the function: 


k{x,y) = 


1 


27r/f y 


exp 



K^(x,y) rf-(x,y)\ 

I 


( 1 ) 


where (v,y) are the offsets from the centre of the Gaussian and 
Ai( and are the standard deviations along (/c, 77 ) = (xsinif/ - 
y cos if/, Xcos if/ - ysinif/), the position-angle-rotated frame of ref¬ 
erence. The FWHMs (a = V 8 In 2 A^ and b = V 8 In 2 Arj) and the 
position angle are usually read from the header as the bmaj, bmin 
and BPA keywords, respectively, but they can also be manually sup¬ 
plied. In addition, for Integral Field Unit (IFU) data, ^^Barolo can 
receive as input an image (or data-cube) of a star, usually observed 
at the same time of the scientific target. The star is then used to 
determine the PSF by fitting it with a 2D Gaussian. 

This smoothing step is the bottleneck of the fitting algorithm, 
since the convolution is a very computational expensive opera¬ 
tion and ^^Barolo needs to perform it for each calculated model. 
To speed-up this step with Fast-Fourier transforms, we used the 
shared-memory parallel OpenMP implementation of the FFTW3 
library. The user is however advised to use data-cubes with suitable 
sampling in order to save computational time. 


2.3 Normalization 

The normalization allows the code to exclude one parameter from 
the fit, namely the surface density E of the gas. We have currently 
implemented two different kinds of normalization: pixel-by-pixel 
and azimuthally averaged. In the former case, the model is nor¬ 
malized such that the column density maps of model and obser¬ 
vations are the same. In other words, we impose that the integral 
of each spatial pixel along the spectral dimension in the model 
is equal to the integral of the correspondent spatial pixel in the 
observations. This type of normalization allows to have a non- 
axisymmetric model in density and avoids that untypical regions, 
like areas with strong and clumpy emission or holes, might affect 
the global fit (see e.g. |Lelli et aL]|2012| ). In the second case, the 
model is instead normalized to the azimuthal-averaged flux in each 
ring. According to our tests, the pixel-by-pixel normalization is of¬ 
ten a more advisable solution, so this is the default in ^^Barolo. 
The azimuthal-averaged normalization is useful to determine the 
inclination angle of the outer rings. The normalization step can be 
turned off, in this case the user can supply a surface density profile 
or leave it free and fit it together with the other parameters. 


^obs. 



Figure 2. Geometrical parameters of the disc model. The galaxy disc in 
the x'y'z' space is projected into an ellipse in the xy plane of the sky. The 
inclination angle i is taken with respect to the line of sight, the position 
angle (p identifies the position of the major axis on the receding half of the 
galaxy and it is taken counterclockwise from the North direction. 


tion algorithm and defines whether a model is suitable or not, is the 
averaged sum of the residuals over each pixel: 


F = I V Ar; w( 6 »;) ( 2 ) 

n 

i=\ 

where n is number of pixels where the residual Ar are evaluated 
and w{6) is a weighting function. To be considered, a pixel must 
either have a non-zero flux in the model or be part of the identified 
emission region of the galaxy. We provide three kinds of residuals: 


Ar = 


(M - Df 

~~W~ 


(3a) 


Ar=\M-D\ 


(3b) 


2.4 Residuals 


\M-D\ 
(M + D) 


(3c) 


The residuals are calculated by comparing the model and the data where M and D are the flux values of the model and the data, re- 

pixel-by-pixel. The number F, which is returned to the minimiza- spectively. The ( [^ residual is a kind of;^^^ without however a con- 
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ventional statistical meaning. When D is a blank pixel, in ( [^ we 
set D equal to the root mean square value (rms) of the cube. The 
( [^ residuals give more weight to regions where the emission is 
faint and diffuse, the ( [^ is intermediate. The weighting function 
in Eq. is w{0) - looser, where 6 is the azimuthal angle (0° for 
the major axis) and m - 0,1,2. With m ^ 0, the residual gives 
prominence to regions close to the major axis, i.e. where most of 
the information on the rotation motions lies. As an option, Eq. ^ 
can be multiplied by a factor (1 + where n\y is the number of 
pixels having the model but not the observation, in order to penal¬ 
ise models that extend farther than the data. This option is useful, 
for instance, to determine the correct inclination of the outer rings. 
As mentioned, this choice should be combined with an azimuthal 
normalization of the surface density. 


2.5 Additional features 

Here, we list and briefly describe some other useful tools available 
inside ^^Barolo: 


- Source det ection. A source finder derived from Duchamp 
( [Whiting 2012| ) is implemented. ^^Barolo can identify all the 
sources in the given data-cube and automatically fit each of them. 

- Masking. In order to obtain a good fit, ^^Barolo needs to build 
a mask to identify the regions that are ascribable to the galaxies in 
the data-cube. The default algorithm takes advantage of the source 
finder and builds a mask directly on the identified emission regions. 
As an alternative, ^^Barolo can smooth the original data-cube by 
a factor n and consider only the regions where the flux is higher 
than m times the rms of the smoothed data-cube. Default values are 
n - 2 and m = 3, but they can be changed by the user. 

- Automatic initial guesses. ^^Barolo can automatically estim¬ 
ate the initial parameters for the fit. The algorithm starts by isol¬ 
ating the galaxy through the source finder. The coordinates of 
the centre are taken as the flux-weighted average positions of the 
source. The systemic velocity is calculated as the midpoint between 
the velocities corresponding to the 20% of the peaks of the global 
line-profile. The position angle is estimated from the velocity held 
as the straight line that maximizes the gradient in velocity along the 
line of sight. The inclination is calculated from the column density 
map: a model map is built from the observed gas density profile, 
smoothed to the same spatial resolution of the data and fitted to the 
observed map. The rotation velocity is calculated as the inclination- 
corrected half-width of the global line profile at 20% of the peak 
flux (W 20 ). Default values for the velocity dispersion and the disc 
thickness are 8 kms“^ and 150 pc, respectively. The algorithm is 
able to determine good initial guesses in most cases. The weakest 
step is the estimate of the inclination, which is also the most critical 
parameter to fit (see also |Begeman|1987] l. These initial estimates al¬ 
low ^^Barolo to be run automatically, but the user can also decide 
to set some parameters and let the code estimate the others. A num¬ 
ber of output plots (e.g. position-velocity diagrams along the major 
and minor axes) are provided to allow the user to check the quality 
of the automatic fits. 


- Regularization of the parameters. Eitting the geometrical 
parameters together with the kinematic ones can lead to unphysical 
discontinuities in the derived rotation curve. The fits of the inclin¬ 
ation and the position angle, in particular, often show unrealistic 
oscillations and numerical scatter. The usual approach for dealing 
with this issue, also in the 2D tilted-ring model, consists in deriv¬ 
ing a first model by fitting all the parameters together, and a second 


model by fixing the geometrical parameters to some functional 
form and fitting only the circular velocity. ^^Barolo can automat¬ 
ically perform this parameter regularization and then proceed to a 
second fitting step, with only rotation velocity and velocity disper¬ 
sion left free. The other parameters are interpolated and fixed either 
to a polynomial of degree m (chosen by the user, default m = 3) or 
a Bezier function. 

- Errors. There is no direct way to calculate the errors on the 
fitted parameters in the 3D approach. We estimate errors via a 
Monte-Carlo method. Once the minimization algorithm converges, 
^^Barolo calculates a number of models by changing the paramet¬ 
ers with random Gaussian draws centred on the minimum of the 
function. This allows the code to oversample the parameters space 
close to the minimum. In this region, the residuals usually have 
the behavior of a quadratic function. Errors for each parameter are 
taken as the range where this quadratic function shows a residuals 
increase of a w percentage with respect to the minimum (default 
value is 5%). Albeit this procedure is not optimal and computation¬ 
ally expensive, it returns errors which are in good agreement with 
those obtained with more standard methods, for example the differ¬ 
ence between the receding and the approaching halves of the disc 
(e.g. |Swaters|1999j l. 


3 TESTING 3DBarolo 

^^Barolo has been already tested on about one hundred galax¬ 
ies, especially local systems observed in the Hi-line, but also lEUs 
(SINEONI and KMOS) data-cubes of high-z systems (Di Teodoro 
et al., in prep). In this section, we show some applications of 
^^Barolo to Hi data-cubes at different spatial resolutions and we 
compare our results with the traditional 2D approach. 


3.1 High-resolution data 

We run ^^Barolo on several high-resolution galaxies from the 
available Hi-survey, like The HI Nearby Galaxy Survey (THINGS, 
[Walter et~^|2008| ), the Very Large Array - ACS Nearby Galaxy 
Survey Treasury (VLA-ANGST, [Ott et al.|[2012[ ) and the Hydro- 
gen A ccretion in LOcal GAlaxies Survey (HALOGAS, |Heald et al.[ 
20ll\ . ^^Barolo has been able to model these high-resolution 
galaxies and the derived rotation curves agree very well with the 
already published 2D rotation curves. 

Here we focus on the automatic modeling of the well-known 
spiral galaxy NGC 5055. This galaxy is a difficult case because 
it has a prominent warp in the outer disc ( [Bosma | [T978] ). We 
used THINGS natural-weighted data-cube at a spatial resolution 
of 10" (about 0.5 kpc, assuming a distance of 10 Mpc) and we 
run ^^Barolo by supplying only initial guesses for the inclina¬ 
tion and the position angle. ^^Barolo automatically identifies the 
galaxy emission in the data-cube, builds the mask, estimates the 
initial guesses for all the other parameters, performs the first fit¬ 
ting step for each radius with all parameters free and the second 
fitting step with only V^ot and CTgas free, after the regularization of 
the other parameters by Bezier functions (inclination and position 
angle) or constants. The comparison between the observations and 
the final best-fit model is shown in Eig. [^through seven represent¬ 
ative channel maps (a standard output of ^^Barolo). In Eig.|^we 
show the resulting values of the fitted parameters in both steps and 
we compare them with those of |de Blok et al.[ ( [T008| ) derived from 
a tilted-ring fit of the velocity field. Note how the regularization of 
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V = 688 km/s v = 626 km/s v = 564 km/s v = 496 km/s v = 434 km/s v = 367 km/s v = 305 km/s 



Figure 3. Comparison between the observations (blue, top panels) and the model (red, bottom panels) for NGC 5055 from the THINGS survey. Here we show 
some representative channel maps taken from a 30" smoothed data-cube. The size of the field is 25' x 25'. The lower contours are at 2.5cr level, negative 
contours are shown in grey. The map corresponding to the systemic velocity of the galaxy is shown in the centre (496 kms“^). Despite the observations have 
an uneven noise distribution, ^^Barolo is able to reproduce the data very well and in particular it traces the prominent warp visible as a deformation of the 
channel maps in the outer regions of the galaxy. 


the geometrical parameters improves the rotation curve, removing 
for instance the unphysical discontinuities around 28 and 47 kpc. 
Our final rotation curve in Figj4jis in very good agreement with 
that derived by |de Blok et al.| { |2008| ) from the same data-cube. The 
main differences arise from kinematic asymmetries between the ap¬ 
proaching and the receding halves of the galaxy (see |Battaglia et al.| 
|2006| ). While the 2D tilted-ring model applied to the entire disc in 
an asymmetric galaxy mainly results in a rotation curve usually av¬ 
eraged between the approaching and the receding side, ^^Barolo 
always finds the model that has the lowest residuals with respect to 
the data. This may result in a rotation velocity that is the average 
between the two halves, or that follows more closely one side rather 
than the other. This is what happens to our rotation curve in Fig.[^ 

The computational time with a regular dual-core laptop for 
running ^^Barolo on this THINGS data-cube, sized 1024x1024x87 
pixels, is about one day. High-resolution observations are also suit¬ 
able to study asymmetries and peculiarities, such as streaming mo¬ 
tions, lopsidedness and extra-planar gas. ^^Barolo always mod¬ 
els galaxies with a single rotating disc and it is not designed to 
handle these peculiarities. In this context, the newest 2D fitting 
codes might be a more desirable choice, since they guarantee lar¬ 
ger possibilities for modeling kinematical and/or geometrical an¬ 
omalies. The TiRiFiC code also provides for large possibilities of 
complex modeling, even for those systems where the 2D approach 
can not be used, like galaxies close to edge-on or having strong 
spiral arms and thick discs (e.g. |Kamphuis et al.||2013| |Schmidt| 
let al.|2014| ). 

3.2 Mid-low resolution data and robustness 

Only a fraction of local galaxies can be observed with a spatial res¬ 
olution comparable to that of the THINGS survey. Most emission¬ 
line observations of galaxies, both from radio-interferometers and 
IFU instruments, currently have less than a dozen resolution ele¬ 
ments throughout the entire disc. In these conditions, the beam 
smearing could heavily affect the derivation of rotation curves with 
2D techniques. ^^Barolo is instead conceived to work with these 
low-resolution data. 

We made a robustness test by using a sample of galaxies from 
the Westerbork Hi survey of Irregular and Spiral galaxies Project 


(WHISP, [^der Hulst, van Albada & Sancisi|2001| l. WHISP com¬ 
prises approximately 350 galaxies in the Local Universe observed 
in Hi with the Westerbork Synthesis Radio Telescope (WSRT). We 
selected 32 galaxies with published reliable rotation curves ( |Swa-| 
|ters|1999| ). The sample only contains dwarf late-type galaxies, with 
Vrot < 100 km s"^ that usually have a poor spatial resolution and a 
relatively low signal-to-noise ratio (S/N). For these galaxies, [Swa-| 
|ters| \\999) derived rotation curves using the following procedure 
that includes a correction for the beam smearing. For each galaxy, 
a first estimate of the rotation curve was determined by interact¬ 
ively fitting the rotation velocity, position angle and inclination as 
a function of radius simultaneously to a set of six position-velocity 
diagrams taken at different angles. This was done by using the In¬ 
spector routine in the GIPSY package, which allows to visually in¬ 
spect different slices and manually tune up the parameters. Centres 
were fixed to the optical values, the systemic velocities were de¬ 
termined from a tilted-ring fit to the velocity field with the centre 
fixed. Then, Swaters constructed a 3D model (Galmod in GIPSY) 
from the estimated parameters and he visually compared it with 
the observations, iteratively adjusting the input rotation curve and 
building a new 3D model until the match between model and ob¬ 
servation was satisfactory. Such a procedure takes a long time for 
each galaxy and it may be subjective, but it was the only way to take 
into account the beam smearing effect and derive reliable rotation 
curves in these low-resolution data. The approach of ^^Barolo is 
analogous, but every step is automatically performed and the best 
model is quantitatively determined. 

We run ^^Barolo in a semi-blind fashion on the sample of 
32 dwarf galaxies. We used 30" smoothed data-cubes and we set a 
ring width of 15", following |Swaters| ( |l999| ). The systemic velocit¬ 
ies and the coordinates of the centre of the galaxies were automat¬ 
ically estimated through the source finding algorithm and fixed to 
those values. In each case, values perfectly compatible with those 
of |Swateri] ( |1999| ) were obtained: the maximum deviations from 
Swaters values are AVsys < 2 kms“^ and Aaq - Ajo < 1 pixel = 
10". Global initial estimates for rotation velocities, velocity disper¬ 
sions and thicknesses of the discs were set to Vrot = 50kms"\ 
cTgas = 8kms"^ and zo = 200 pc for all galaxies, while the ini¬ 
tial inclinations and the position angles were taken from Tab. Al, 
Chap. 4 of |Swater^ ( |1999| ). We used the | M - D | residuals with a 
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Figure 4. Tilted-ring parameters derived with ^^Barolo using the THINGS data-cube of the warped galaxy NGC 5055. Left panels: from the top, the rotation 
curve, the inclination and the position angle for each ring. Right panels: from the top, the velocity dispersion, the systemic velocity and the coordinates of the 
centre for each ring. The grey-empty dots represent the first fit with all parameters kept free. The red dots or lines represent the second fit of Vrot and cTgas after 
the regularization of the other parameters (Bezier interpolation for i and constant value for V^ys, xq and yo)- The cyan-shadowed region and the cyan-dashed 


lines are the rotation curve (with errors) and the geometrical angles derived by 

de Blok et al.|(|2008) with a 2D fit. For their final rotation curve, de Blok et al. 

used the coordinates of the centre (xo,yo) = (517,512) from Trachternach et al. 

J2008| and a Fsys = 496.8 km s“^, in perfect agreement with our values within 


the errors. 


cos(^) weighting function. Masks were built by smoothing the data 
by a factor 2 and considering only those regions with flux > 3rms, 
being rms the root mean square (noise) of the smoothed data-cube. 
^^Barolo performed a first step by fitting V^ot , cri and 0 and 
a second step by fitting only V^ot and cTgas and fixing i and 0 to a 
2 nd polynomial function. The execution time is less than a 

minute per galaxy on a regular laptop. 

A comparison between Swaters’ rotation curves and ours is 
shown in Fig. Position-velocity diagrams (output of ^^Barolo) 
along major and minor axes are shown in Appendix A (available 
as online material). In general, the agreement is very good. Most 
differences in the rotation velocities can be attributed to asymmet¬ 
ries in the kinematics between the receding and approaching halves 
of the galaxies. It is interesting to notice that in some galaxies our 
rotation curves rise more steeply than Swaters’ in the inner regions 
(e.g. UGC5272, UGC 6446, UGC 9211). Since the main effect of 
the beam smearing is to reduce the velocity gradients in the rising 
part of the rotation curve, it is possible that the correction manually 
made by |Swater'^ ( |1999| ) was not sufficient in those cases. 

Out of 32 data-cubes, ^^Barolo failed in determining accept¬ 
able models for 4 galaxies, either not converging or deriving wrong 
kinematics. Two cases (UGC 3966, UGC 8837) can be attributed to 


a wrong fit of the inclination in the first step and fixing it to the ini¬ 
tial value led to a good model (see Fig. [^. In the last two cases, 
UGC 7690 and UGC 8490, we could not make the code work¬ 
ing manually neither. The first galaxy is very faint and the fit is 
hampered by the noise. The galaxy UGC 8490 has a huge warp in 
inclination and the algorithm tries to reproduce it by varying the ve¬ 
locity rather than the inclination angle. The final model looks good, 
but the rotation curve is probably unphysical. In addition, from the 
inspection of the position-velocity diagrams along the minor axis, it 
turned out that the automatic estimate of the centres was slightly in¬ 
accurate in five galaxies. Putting the optical centres manually led to 
better models, even though the rotation curves did not change sig¬ 
nificantly. Overall, 94% of the galaxies were accurately modeled 
by ^^Barolo. 

From this test emerges that ^^Barolo is able to derive reli¬ 
able kinematics in low-resolution and noisy data-cubes. We remind 
that we run ^^Barolo in a almost blind execution, since the only 
information we supplied to the code were the initial guesses for 
the inclination and the position angles. Of these, the inclination is 
especially critical as it may be often unsuccessfully estimated by 
the code. Improvements in the initial parameters estimate algorithm 
will be considered in the next releases. We stress however that the 
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Figure 5. Rotat ion curves of the 30 dwarf late-type galaxies selected from the WHISP sample. The gray-shadowed regions represent the rotation curves from 
|Swaters||l999t within the errors, the red dots are the rotation curves derived with ^^Barolo. Since Swaters’ errors are symmetric, his points would lie in 
the centre of the grey band at the same radii that our points. In UGC 7399, our rotation curve stops earlier because in our data-cube there was no signihcant 
emission beyond 175 arcsec. 
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careful inspection of the outputs (position-velocity diagrams and 
model cubes) does clearly single out cases where the fit was not 
successful. Our test on the WHISP data reveals that the success 
rate is very high. This, combined with the very low computational 
time needed to fit these low-resolution galaxies (< 1 minute) are 
key features for the application of ^^Barolo to the upcoming large 
Hi surveys. Indeed, already planned Hi survey, such as WALLABY 
and DINGO with ASKAP {[Johnston et al.|2008|l, LADUMA with 
MeerKAT ( [Booth et al.|2009) and WNSHS with WSRT/APERTIF 
jVerheijen et al.|2008) , are expected to observe thousands of galax- 
ies with spatial resolution comparable to the WHISP galaxies. 


3.3 Very low-resolution data 


In this section we show the effect of the beam smearing on the 
derivation of rotation curves from 2D and 3D analysis, going down 
to extremely low spatial resolution. 

We run both ^^Barolo and Rotcur on a small sample of 
nearby galaxies observed in Hi with single-dish telescopes. We 
selected 4 galaxies (NGC 2403, NGC 2903, NGC 3198 and 
NGC 5055) for which we have both very high resolution Hi data 
( jFraternali et al.j[20021 and THINGS) and single-dish observa¬ 
tions. Effelsberg data for NGC 3198 and NGC 5055 were kindly 
provided by B. Winkel (jWinkel, Floer, & Kraus||2012| |Winkel,| 
[Kraus, & Bach|2012> . NGC 2403 <de Blok et al.|20T^ and NGC 
2903 (Pisano et al., in prep.) where observed with the Green Bank 
Telescope (GBT) and kindly supplied by D. J. Pisano. We com¬ 
pared high-resolution rotation velocities and velocity dispersions 
obtained with ^^Barolo with the low-resolution ones derived both 
with 2D and 3D approaches. Typical spatial resolution is 6" for the 
high-resolution data and 650" (Effelsberg) or 525" (GBT) for the 
low-resolution data, which means that these galaxies are barely re¬ 
solved. High-resolution rotation curves and dispersions were ob¬ 
tained with ^^Barolo from the natural-weighted THINGS data- 
cubes as described in Section [TT] Low-resolution velocity and dis¬ 
persion fields were derived as and 2^^ moment maps, respect¬ 
ively (see e.g. NGC 2403, Fig. |^. In the 2D approach, rotation 
curves were derived with Rotcur fitting only V^ot and fixing the 
other parameters to the high-resolution values, 2D velocity disper¬ 
sion profiles were obtained by taking the average value along the 
rings on the dispersion fields. No correction for the beam smear¬ 
ing was performed. ^^Barolo was run by fitting only Vrot and CTgas, 
except for NGC 5055, where we kept free also the position angle 
to trace the outer warp. For NGC 2403, we used only the reced¬ 
ing half of the disc, since the approaching half is contaminated by 
Hi emission from the Milky Way. In Fig. we show the result¬ 
ing rotational velocities, velocity dispersions and the comparison 
between the models and the data through the position-velocity dia¬ 
grams along the major and the minor axes. As expected, a 2D ap¬ 
proach is not suitable for data at these resolutions, since the beam 
smearing significantly affects the derivation of the 2D maps from 
the data-cubes. Beam smearing flattens the gradients in the velo¬ 
city profiles and turns rotation velocity into apparent high velocity 
dispersion, as it clearly appears from the maps in Fig.[^ Such a de¬ 
generacy between Vrot and cTgas is broken in ^^Barolo, because the 
beam s meari ng effect is taken into account in the convolution step 
(Section ^i. Unlike the 2D approach, ^^Barolo recovers correct 
rotation velocities for all these galaxies in every location of the disc. 
The differences between the low-resolution and the high-resolution 
velocities are within the errors in all cases. Even more remarkably, 
^^Barolo returns low values for the intrinsic velocity dispersions 
of the gas that are fully comparable with the correct values (bottom 
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Figure 6. Velocity fields (L^ moment) and velocity dispersion fields (2^^^ 
moment) for NGC 2403 derived from high (/g/r, jFraternali et al.|2002) and 
low resolution (right, jde Blok et al. pout data-cubes. Images are on the 
same spatial and velocity scale. Note the dramatic effect of the beam smear¬ 
ing: the velocity field at high resolution, showing the typical traits of a fiat 
rotation curve, turns into a nearly solid-body pattern at low resolution, es¬ 
pecially in the inner parts. Velocity dispersions increase by a factor 3-4 
throughout the whole disc. 


panels of the upper plots). An inspection of the position-velocity 
diagrams along the minor axis (bottom panels of the lower plots) 
should give an idea of how the line broadening in these data is fully 
dominated by instrumental effects. 

These results show that ^^Barolo is a powerful tool to study 
the kinematics of galaxies even in very low-resolution data (2-3 
resolution elements across the whole disc), where the standard 2D 
approach fatally fails. Above all, ^^Barolo is almost always able 
to describe the correct shape of the rotation curves and it becomes 
extremely robust when its outputs are visually checked by the user. 
A particularly interesting application concerns the Integral Field 
Spectrographs (IFS) data. Instruments like SINFONI ( jEisenhau^ 
et al. 2003j ), KMOS ( [Sharpies et al.|200^ and MUSE ( [Bacon et~ 
20101 on the Very Large Telescope (VLT) can observe line emis¬ 
sion, such as Hex, N and O forbidden lines, in galaxies up to red- 
shift about 2.5, with a resolution similar to the observations showed 
in this section. Running ^^Barolo on those data-cubes is a very 
challenging task because, in addition to the poor spatial resolution, 
these observations have a low spectral resolution (channel widths 
of 30-40 kms"^ compared to the few kms"^ of Hi data) and very 
low S/N. Nevertheless, a 3D approach should be heartily recom¬ 
mended in order to take advantage of the full information available 
in these data-cubes. 
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Figure 7. Upper panels: velocities (top) and velocity dispersions {bottom) for NGC 2403, NGC 2903, NGC 3198 and NGC 5055. Blue dots were derived from 
THINGS high-resolution (6") data-cubes with ^^Barolo. Cyan regions represent the errors. Red triangles are the fit with ^^Barolo in single-dish data, the 
green open diamonds are the results obtained from the 2D maps (Rotcur on velocity fields and Ellint on dispersion fields). Lower panels: the correspondent 
position-velocity diagrams along major {top) and minor {bottom) axes. Data are represented in grey, models in red, rotation curves as blue square dot. 


4 TESTING ^i^BAROLO’S LIMITS 

The goodness of a fit is mainly determined by the combination of 
four factors: the inclination of the galaxy disc, the spatial resolu¬ 
tion, the spectral resolution and the S/N. We used mock galaxies to 
find some significant thresholds to these quantities and test under 
what conditions ^^Barolo may have problems in deriving a reliable 
kinematic model. 

We built initial artificial data-cubes with a resolution of 12" 
(FWHM), a channel width of 7.5 kms“^ and no noise. The pixel 
size is 3", thus the beam covers an area of about 18 pixels. The 
galaxy models have a maximum radius of 400", i.e. about 33 spa¬ 
tial resolution elements per side. This configuration could be a typ¬ 
ical observation of a nearby galaxy with a modern interferometer, 
like the JVLA. We set global parameters for aq, yo» ^sys and zo, 
an exponential gas density profile and a constant dispersion field 
of 10 kms"^ The initial inclination is i - 60°. We assumed ro¬ 
tation curves shaped as Vrot(^) = 2/:7r Vq arctan (/?//?o). being Vq 


the asymptotic circular velocity and Rq the turnover radius, i.e. the 
transitional point between the rising and fiat part of the rotation 
curve. We made a model with a steeply-rising plus fiat rotation 
curve (Vo = 250 kms"\ Rq - 20") and a model with slowly-rising 
solid body-like rotation curve (Vo = 150 kms"\ Rq = 150"). We 
progressively degraded the artificial data-cubes and we try to get 
back the input rotation curves using ^^Barolo. 

We first reduced the spatial resolution by smoothing the ini¬ 
tial data-cubes down to 200" (preserving the number of pixels per 
beam), that is 2 resolution elements per side for the model galaxies. 
With a large number of velocity channels and no noise, ^^Barolo 
is able to recover the correct rotation curve even at the lowest res¬ 
olution. This could be the case of a single-dish observation as those 
in Section [33] For comparison, we also derived the rotation curves 
running Rotcur on the Gaussian velocity fields, fitting only the ro¬ 
tation velocity. Both approaches recover an almost perfect rotation 
curve when the galaxy has more than 10 resolution elements, but, 
below this limit, Rotcur increasingly underestimates the rotation 
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Figure 8. Rotation velocity and velocity dispersions derived with ^^Barolo 
from models with flat (top) and slowly-rising (bottom) rotation curves. Here 
we show the lowest resolution models, i.e. 100" (left) and 200" (right), and 
three inclinations (45°, 60° and 75°). Models are the black thick lines. At 
these resolutions, we run ^^Barolo with just V^t and crg^s free. 


velocity, especially in the inner regions of the model with flat rota¬ 
tion curve, whereas ^^Barolo can successfully determine it. With 
the 2D approach the relative errors with respect to the actual rota¬ 
tion curve are up to 60% in the lowest resolution and in the inner 
parts, whereas with the 3D approach the errors are confined to a few 
percent at all resolutions and over the entire disc. The derived rota¬ 
tion curves obtained at 100" and 200" (4 and 2 resolution elements 
respectively) are shown in Fig.[^(red circles). 

Next, we studied the effect of the inclination angle on 
^^Barolo’s accuracy, in particular on the fit of Vrot and cTgas. We let 
the inclination of the model galaxies varying from nearly face-on to 
nearly edge-on and fitted these models by fixing the inclination. In 
Fig. we show the recovered rotation curves and velocity disper¬ 
sions for the lowest resolution models, namely two and four resol¬ 
ution elements per side, and for three significant inclination angles, 
i.e. 45°, 60° and 75°. The rotation velocity is well recovered at any 
inclination both in models with fiat and solid-body rotation curve, 
although in the flat model for / > 75° the inner points of the rota¬ 
tion curves start to be underestimated. The velocity dispersion of 
the inner point can be overestimated in some cases by a factor up 
to about 2 but with large error bars, the other points are recovered 
within a few km s ~^. 

Overall, if the inclination is known, ^^Barolo can derive the 
correct rotation curve and can disentangle between rotation and ve¬ 
locity dispersion at almost every inclination and even for data at 
very low spatial resolution where the 2D approach can not be used. 
When the inclination is not known, the code can fit it together with 
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Figure 9. Relative percentage errors for the fitted Vrot (red circles) and cTgas 
(cyan diamonds) as a function of the average S/N at fixed i = 60°. We show 
the results for the flat (top) and solid-body model (bottom) at the lowest 
resolutions, i.e. 100" (4 resolution elements, left) and 200" (2 resolution 
elements, right). The code becomes unreliable below a S/N < 2-2.5. 


the kinematical parameters, but this requires some care. The inclin¬ 
ation is the thorniest parameter to deal with, and running ^^Barolo 
with a completely unknown inclination is not advisable since the 
code can easily converge to a local minimum close to the initial in¬ 
clination. The initial guess for the inclination is therefore essential 
for the goodness of the fit and it can be either supplied by the user 
(e.g. optical values) or estimated by ^^Barolo. Our tests show that 
the algorithm for the initial guesses (see Section [O] ) returns good 
estimates (within a few degrees) of the global inclination in most 
cases regardless of the spatial resolution when 45° < i < 75°. Thus, 
in general for i > 45° one should always be able to obtain accept¬ 
able kinematical fits (note that above 75°, errors in the inclination 
have little impact in the rotation velocity). However, for i < 45° ro¬ 
tation curves may become progressively more uncertain due to the 
smaller rotational component along the line of sight and the larger 
impact of inclination errors. This is a problem for any fitting al¬ 
gorithm. Finally, we note that the inclination is degenerate with the 
disc thickness and we expect this effect to be important especially 
for dwarf galaxies. In the future, we will consider to include a self- 
consistent treatment (assuming the hydrostatic equilibrium of the 
gas) for the disc thickness in ^^Barolo. 

We also tested the effect of the noise on the derivation of ro¬ 
tation velocities and velocity dispersions. We focused on i - 60°, 
used a constant surface density profile for the galaxy models and 
added a progressively higher Gaussian noise to the artificial data- 
cubes with two and four resolution elements. Relatively low S/N 
and less than 10 resolution elements characterize ALMA data of 
high-z galaxies and are expected for most data of the future Hi sur¬ 
veys. We calculated the relatives percentage errors of the fitted V^ot 
and CTgas with respect to the true values. In Fig.we show the aver- 
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age errors over the entire disc as a function of the average S/N. The 
performance of ^^Barolo remains consistently high whenever the 
signal to noise of the average emission in the ring is S/N > 2. Be¬ 
low this limit, ^^Barolo may consider portions of the background 
as galaxy emission, leading to a wrong fit of the circular velocity 
and velocity dispersion, which can be either overestimated or un¬ 
derestimated without systematic changes in slope or shape. Further 
tests in the inclination range 45° < i < 75° did not significantly 
change such a threshold and the trends showed in Fig.|^ We stress 
that these tests were run with the default masking options and that 
a manual fine-tuning of the mask might lead to good results even at 
slightly lower signal-to-noise ratios. 

Finally, we reduced the spectral resolution in the noisy artifi¬ 
cial data-cubes. Low spectral resolutions characterize, for instance, 
IFUs data (like KMOS or MUSE), which typically have channel 
widths of 30-40 km s"^ for each emission line, i.e. observed galax¬ 
ies can span less than ten channels. Increasing the channel width 
dramatically lowers the number of data points that ^^Barolo can 
use to constrain the best model. From our test emerges that the 
number of channels that guarantees a good fit at low S/N varies 
between 8 and 12, depending on the spatial resolution. 

In conclusion, ^^Barolo can work even with observations at 
very low spatial/spectral resolution and low S/N and in a wide 
range of galaxy inclinations. These factors influence together the 
goodness of the model. In an extreme case of a galaxy with just a 
couple of resolution elements per side, it would be advisable that 
the source is detected at a S/N > 3 and over a dozen channels or 
more. 


5 CONCLUSIONS AND FUTURE PROSPECTS 


In this paper, we presented ^^Barolo, a new software for the fit¬ 
ting of 3D tilted-ring models to emission-line observations, and we 
showed several examples of applications, from high-resolution to 
very-low resolution data-cubes. The main purpose which led us 
to develop ^^Barolo is to have a tool to derive the kinematics of 
galaxies using a simple rotating disc model, feasible to be used at 
a wide range of resolutions without being affected by instrumental 
effects. At high resolution, ^^Barolo works as well as the 2D tradi¬ 
tional approach, but it is more computational expensive and it does 
not yet allow to study the peculiarities of galaxies, such as non¬ 
circular motions. ^^Barolo reaches its best performance in deriving 
reliable rotation curves from low-resolution data, down to barely 
resolved galaxies, where the 2D approach totally fails because of 
the beam smearing. Moreover, the robustness and the possibility 
of identifying sources and estimating the initial conditions for the 
fit make ^^Barolo a tool that can be automatically run on large 
amounts of data. 

^^Barolo can operate on all emission-line data. In particular, 
suitable observations are the Hi-line, CO and C^ lines from sub¬ 
millimeter interferometers, such as the Plateau de Bure Interfero¬ 
meter (PdBI) or ALMA, and optical/IR recombination lines from 
the last generation of IFUs. These instruments are providing the 
extraordinary opportunity to study the evolution of the kinematics 
and the dynamics of galaxies throughout the Hubble time. Several 
recent studies have used these data to investigate how the earliest 
disc galaxies differ from local galaxies through different analysis, 
like the V/cr ratio, which measures the dynamical support given by 
rotation (e.g. [Forster Schreiber et al.|2()0^|Mancini et al.|201 1| ), the 
Toomre Q p arameter (e.g.jGenzel et al.|2011| t or the turbulence in 
the disc (e.g. Green et al. 201 Oj ) . Running B arolo on these obser¬ 


vations will give the opportunity to derive reliable rotation curves 
at high-z and break the degeneracy between rotation velocity and 
velocity dispersion down to very low spatial resolutions. 

Overall, the kinematics is the starting point to investigate sev¬ 
eral properties of disc galaxies, including but not limited to the 
dynamics, the matter distribution and the evolution of scaling re¬ 
lations (e.g. [Epinat et al.|2009[ [Miller et al.|2012| l. The ability to 
quickly derive rotation curves will be more and more important in 
the future given the increasing number of emission-line observa¬ 
tions available from various astronomical facilities. 
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APPENDIX A: P-V DIAGRAMS OF WHISP GALAXIES 

In this appendix we show the comparison between our best mod¬ 
els and the observations for 30 dwarf galaxies selected from the 
WHISP sample (see Section 3.2). Ror each galaxy we show slices 
along major and minor axes (outputs of ^^Barolo). Data are in 
grey and cyan (negative contours) and the model in red. Dark-blue 
square are the derived rotation curves. 
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